module mod_elgmc
use mod_elglike
use mod_elgdata
implicit none

contains
  
  subroutine elgmc_mcmc(in_params, steps, filename) 
    type(t_params), intent(inout) :: in_params
    type(t_params)  :: params

    type(t_params) :: proposed
    integer, intent(in) :: steps
    type(t_elglike) :: elglike
    double precision :: L1, L2, accept, pr, rand
    integer :: file, i, lmax
    character(len=128), intent(in) :: filename
    

    file = 3

    call params_copy(in_params, params)

    lmax = global_lmax
    call elglike_initialize(elglike, lmax)

    write(*,*) "Calculating first step.."
 !   call elglike_calculate_likelihood_cheat(elglike, params)

   call elglike_calculate_likelihood_master(elglike, params)
    L1 = elglike%likelihood

    open(unit=file,file=filename,action='write',status='replace')
    i = 0
    call params_savetext(params,file,i, filename)
    
    close(file)


    accept = 0
    do i=1, steps
       
       if (mod(i,250) ==0 ) then
          params = in_params
       end if
          
       write(*,*) "Step",i, "..."
       call params_propose(params, proposed) 
       call params_print(proposed)
       !write(*,*) "steplength", params

       call elglike_calculate_likelihood_master(elglike, proposed)
    
       !call elglike_calculate_likelihood_cheat(elglike, proposed)
    
       L2 = elglike%likelihood
       
       Pr = l2/l1
       
       !write(*,*) pr, l1, l2
       if (Pr>=1.0) then
          call params_copy(proposed, params)
          accept = accept + 1
          l1 = l2
       else 
          call random_number(rand)
          !write(4,*) i, rand
          if (rand<Pr) then
             call params_copy(proposed, params)
             accept = accept + 1
             l1 = l2
             !write(*,*) "MCMC step accepted"
          else
             write(*,*) 'SAMPLE REJECTED'
          end if
       end if
       call params_savetext(params,file,i, filename)
       write(*,*) 'acceptance ratio = ', real(accept) / real(i)
    end do
!    close(4)
    
  end subroutine elgmc_mcmc


end module mod_elgmc
